// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "sparse.h"

template<typename Scalar>
void
initSPD(double density, Matrix<Scalar, Dynamic, Dynamic>& refMat, SparseMatrix<Scalar>& sparseMat)
{
	Matrix<Scalar, Dynamic, Dynamic> aux(refMat.rows(), refMat.cols());
	initSparse(density, refMat, sparseMat);
	refMat = refMat * refMat.adjoint();
	for (int k = 0; k < 2; ++k) {
		initSparse(density, aux, sparseMat, ForceNonZeroDiag);
		refMat += aux * aux.adjoint();
	}
	sparseMat.setZero();
	for (int j = 0; j < sparseMat.cols(); ++j)
		for (int i = j; i < sparseMat.rows(); ++i)
			if (refMat(i, j) != Scalar(0))
				sparseMat.insert(i, j) = refMat(i, j);
	sparseMat.finalize();
}

template<typename Scalar>
void
sparse_solvers(int rows, int cols)
{
	double density = (std::max)(8. / (rows * cols), 0.01);
	typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
	typedef Matrix<Scalar, Dynamic, 1> DenseVector;
	// Scalar eps = 1e-6;

	DenseVector vec1 = DenseVector::Random(rows);

	std::vector<Vector2i> zeroCoords;
	std::vector<Vector2i> nonzeroCoords;

	// test triangular solver
	{
		DenseVector vec2 = vec1, vec3 = vec1;
		SparseMatrix<Scalar> m2(rows, cols);
		DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);

		// lower - dense
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
		VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
						 m2.template triangularView<Lower>().solve(vec3));

		// upper - dense
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
		VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
						 m2.template triangularView<Upper>().solve(vec3));
		VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
						 m2.conjugate().template triangularView<Upper>().solve(vec3));
		{
			SparseMatrix<Scalar> cm2(m2);
			// Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr
			MappedSparseMatrix<Scalar> mm2(
				rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(), cm2.valuePtr());
			VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
							 mm2.conjugate().template triangularView<Upper>().solve(vec3));
		}

		// lower - transpose
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
		VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
						 m2.transpose().template triangularView<Upper>().solve(vec3));

		// upper - transpose
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
		VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
						 m2.transpose().template triangularView<Lower>().solve(vec3));

		SparseMatrix<Scalar> matB(rows, rows);
		DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);

		// lower - sparse
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular);
		initSparse<Scalar>(density, refMatB, matB);
		refMat2.template triangularView<Lower>().solveInPlace(refMatB);
		m2.template triangularView<Lower>().solveInPlace(matB);
		VERIFY_IS_APPROX(matB.toDense(), refMatB);

		// upper - sparse
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular);
		initSparse<Scalar>(density, refMatB, matB);
		refMat2.template triangularView<Upper>().solveInPlace(refMatB);
		m2.template triangularView<Upper>().solveInPlace(matB);
		VERIFY_IS_APPROX(matB, refMatB);

		// test deprecated API
		initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
		VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
						 m2.template triangularView<Lower>().solve(vec3));

		// test empty triangular matrix
		{
			m2.resize(0, 0);
			refMatB.resize(0, refMatB.cols());
			DenseMatrix res = m2.template triangularView<Lower>().solve(refMatB);
			VERIFY_IS_EQUAL(res.rows(), 0);
			VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
			res = refMatB;
			m2.template triangularView<Lower>().solveInPlace(res);
			VERIFY_IS_EQUAL(res.rows(), 0);
			VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
		}
	}
}

EIGEN_DECLARE_TEST(sparse_solvers)
{
	for (int i = 0; i < g_repeat; i++) {
		CALL_SUBTEST_1(sparse_solvers<double>(8, 8));
		int s = internal::random<int>(1, 300);
		CALL_SUBTEST_2(sparse_solvers<std::complex<double>>(s, s));
		CALL_SUBTEST_1(sparse_solvers<double>(s, s));
	}
}
